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A multidimensional numerical simulation framework aimed at investigating the process of plasma detach- 
ment from a magnetic nozzle is introduced. An existing numerical code based on a magnetohydrodynamic for- 
mulation of the plasma flow equations that accounts for various dispersive and dissipative processes in plasmas 
was significantly enhanced to allow for the modeling of axisymmetric domains containing three-dimensional 
momentum and magnetic flux vectors. A separate magnetostatic solver was used to simulate the applied mag- 
netic field topologies found in various nozzle experiments. Numerical results from a magnetic diffusion test 
problem in which all three components of the magnetic field were present exhibit excellent quantitative agree- 
ment with the analytical solution, and the lack of numerical instabilities due to fluctuations in the value of V ■ B 
indicate that the conservative MHD framework with dissipative effects is well-suited for multi-dimensional 
analysis of magnetic nozzles. Further studies will focus on modeling literature experiments both for the pur- 
pose of code validation and to extract physical insight regarding the mechanisms driving detachment 


Nomenclature 


B 

magnetic induction, T 

P 

pressure, Pa 

Bm 

Maxwell stress tensor 

q 

dissipative power per area, W/m 1 2 

E 

electric field strength, V/m 

T 

temperature, K 

E' 

electric field in the plasma frame, V/m 

u 

fluid velocity, m/s 

£ 

energy density of the plasma, J/m 3 

V 

resistivity, O.m 

k 

thermal conductivity, W/m.K 

po 

permeability of free space, H/m 

j 

current density, A/m 2 

p 

mass density, kg/m 3 


I. Introduction 

M AGNETIC nozzles are used in many experiments to control - confine, cool, accelerate, and direct - plasma 
flows using magnetic fields in lieu of material boundaries. Studies of plasma detachment have been motivated 
by its role in many astrophysical phenomena 1-3 and by the various terrestrial applications that use magnetic nozzles 
to exhaust power and particles. These applications include plasma flows in fusion systems 4 where understanding 
detachment is important to controlling the thermal loads on components, rapid plasma cooling for laser emission 
where understanding recombination during the expansion is crucial to the generation of a population inversion, and 
directional control of a plasma jet in etching processes. 5 We focus on the role an applied magnetic field and magnetic 
nozzle can play in facilitating the transfer of random thermal energy into directed kinetic energy in a plasma thruster. 

Since magnetic fields are solenoidal, the design of a magnetic nozzle for a plasma thruster involves conflicting 
goals of confinement (leading to energy conversion and plasma acceleration along the magnetic field) and detachment 
(allowing for plasma motion transverse to the magnetic field leading to momentum transfer and thrust production). For 
a thruster it is important to understand the process of plasma detachment from an applied magnetic field both in terms 
of developing a capability to accurately model the process and for the development of a basic physical understanding 
of the process, with the end goal being a better understanding of thrusters that employ magnetic nozzles. 
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While there have been many experimental s " 16 and analytical 17-21 efforts to investigate detachment mechanisms, 
multidimensional numerical investigations have been sparse. 22 Plasma properties in real devices can vary significantly 
in space and time. Such variations, as well as the complicated magnetic field topologies found in many nozzles, 
make analytical treatments of detachment difficult if not wholly intractable. In addition many atomic and plasma 
kinetic phenomena, such as the values of various transport coefficients or the non-ideal relationship between the 
pressure, density, and temperature as well as their relationship with the ratio of specific heats in the plasma, cannot be 
easily measured with common experimental diagnostics but can greatly affect the detachment process. Consequently, 
numerical simulations that incorporate realistic models of the various plasma processes can play a unique role in 
developing an understanding of the mechanisms that cause detachment and learning how these processes scale with 
various plasma and applied field properties. 

In this paper, we discuss the recent development of multidimensional numerical tools that can be used to investigate 
plasma detachment mechanisms. The numerical code is based on a conservation formulation of the magnetohydro- 
dynamic (MHD) equations, and contains models for various classical and anomalous transport phenomena and a real 
(non-ideal) equation of state. The code was introduced in Ref. [23] and has been previously used to simulate self-field 
magnetoplasmadynamic thrusters (MPDT) 24,25 and a lithium Lorentz force accelerator. 26 The capabilities of the code 
are significantly enhanced over those previous studies, resulting in the capability to model plasma flows in magnetic 
nozzles. In addition, the code could also be used to simulate applied-field MPDTs, which incorporate the processes 
found in an MPDT with the magnetic nozzle field topology. We present results from a magnetic diffusion test case 
involving multi-dimensional magnetic fields. These results suggest that certain difficulties encountered when applying 
other MHD simulation codes to the magnetic nozzle problem may not arise in our numerical studies. In future studies 
the code will be employed to model various plasma flows that show potential propulsive benefits. The modeling results 
will be compared with available literature data to explain various experimental observations and gain insight into the 
various physical processes associated with plasma flow in a magnetic nozzle and plasma detachment. 

The outline for the rest of this paper is as follows. In Sect, n we review various proposed mechanisms that could 
lead to plasma detachment. The physical model, governing equations, and numerical formulation for the nozzle sim- 
ulations are described in Sect. HI. We present the results of a plasma simulation test case involving multi-dimensional 
magnetic fields in Sect. IV, comparing the analytical solution to a magnetic diffusion problem with corresponding 
numerical simulation results. 


II. Review of Plasma Detachment Theories 

The problem of plasma detachment from a magnetic field has been addressed in various branches of plasma 
physics. A survey of the literature suggests several mechanisms that could lead to detachment. These are listed and 
discussed in turn below. 

A. Resistive Detachment 

A plasma that is sufficiently resistive can support collisionally driven plasma transport transverse to the magnetic 
field. 20 Collisions that accomplish this arise either from classical Coulomb interactions between electrons and ions or 
through anomalous transport driven by interactions between electrons and electromagnetic waves (i.e. plasma turbu- 
lence 27 ). Hooper 20 analyzed detachment for such a case and concluded that the resistively-driven change in magnetic 
flux due to cross-field transport resulted in a violation of the “frozen-in” flow condition leading to a loss of magnetic 
flux over time given by, 

( 1 ) 

where ip is the magnetic flux, e 0 is the permittivity of free space, rj is the plasma resistivity, w p is the plasma frequency, 
u is the plasma velocity, and ^ co u is the effective electron collision frequency. This.relationship can be used to estimate 
the minimum required collision frequency for plasma detachment. 

B. Kinetic Detachment 

Studies of solar wind in the late 1950s led to the proposition that plasma detaches from a magnetic field when the 
flow velocity exceeds the local Alfv6n speed. 28 This corresponds to the case where the kinetic energy of the plasma is 
greater than the magnetic field energy. 29 As the plasma expands in the downstream portion of a magnetic nozzle, the 
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ratio of plasma kinetic energy to the energy in the magnetic field continues to increase. This ratio can also be written 
as the square of the ratio of plasma velocity to the Alfv&n speed v a , 


P • 


pu 2 / 2 / u 

W( 2/*.) ~ VB/ViCP 


-Y 

VAj 


= Ml. 


( 2 ) 


When (3 exceeds 1, it is theorized that plasma will possess sufficient kinetic energy to permit detachment from the 
magnetic field. This theory was recently applied to plasma thrusters by Arefiev and Breizman. 21 Using MHD theory, 
they concluded that plasma would detach from the nozzle, even in the absence of collisions, with an detachment 
efficiency (ratio of the outgoing to incoming axial momentum fluxes) rjd, given by, 



where 0 o is the nozzle divergence angle. 


( 3 ) 


C. Recombination Detachment 

Applied magnetic field confinement of particles in an expanding jet is lost when the recombination rate is such that 
ionized particles are quickly converted into neutrals. 30 However, based on expansion and cooling rates in typical 
thruster flows, it appears unlikely that the recombination rate would be high enough to support this detachment mech- 
anism. It has been suggested 30 that injecting a cold gas of heavy molecules (such as Kr) into the expanding plasma 
flow for the purposes of cooling the electrons and enhancing detachment could lead to enhanced recombination-driven 
detachment. 


D. Non- Adiabatic Detachment 

If a magnetic field has a strong spatial gradient, the condition of adiabatic invariance of a charged particle’s magnetic 
moment may not hold. Specifically, if the ratio of a particle’s gyro radius to the characteristic length scale of change 
in the magnetic field, r g /(B/VB), is greater than unity then charged particles are not confined to follow magnetic 
field lines through invariance of their magnetic moments and detachment can occur. In the types of plasmas found 
in magnetic nozzles, this ratio could be greater than unity for ions but is typically much less than unity for electrons. 
With electrons confined to magnetic field lines, it is unlikely that this detachment mechanism plays a major role in 
propulsive magnetic nozzles since the plasma expelled in these systems is typically quasineutral. 20 

£. Electron Inertia Detachment 

It was proposed that because electrons possess finite mass, they can inhibit the flow of azimuthal currents in the nozzle 
leading to collisionless detachment. 17, 19 A theoretical treatment showed that plasma detachment could occur when the 
non-dimensional parameter, 

G ~ w e uii . ( 4 ) 

was below a critical value. Here w e and Wj are the gyro frequencies of the electrons and ions, and r a and u 0 are 
the radius and velocity of the plasma jet at the throat. This effect, however, is quite minimal and could be taken to 
represent a lower bound on the amount of detachment possible in a given experiment. 

III. Physical Models for Magnetic Nozzle Simulations 

Among the abovementioned detachment theories, the resistive and kinetic mechanisms are most amenable to in- 
vestigation using models based on the continuum treatment of plasma flows. An analysis of detachment driven by 
recombination, non-adiabaticity, or electron inertia requires either a particle-in-cell (PIC) or hybrid treatment of the 
plasma. The numerical framework we describe below is based on the MHD (continuum) formulation of the plasma 
flow governing equations, restricting our studies of detachment to the resistive and kinetic detachment mechanisms. 
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A. MHD Equations 


The governing equations for an MHD flow problem represent the conservation relations for mass, momentum, mag- 
netic flux, and energy. These can be compactly written in the form, 


dt 
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. (£ + p) u - Bm • u . 
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(5) 


The right hand side of the equation, D, contains dissipative effects, which are also written as the divergence of fluxes. 26 

The compact vector-tensor equations are expanded into their individual vector components and rewritten in a 
finite-volume form so that they can be solved by a conservative MHD code. The expanded equation set, assuming no 
azimuthal variation ( d/d& = 0), is: 
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where 


q r = ( E' z B e - E' g B z ) jp Q + ( k rr dT e /dr ) + (k rz dT e /dz ) , 
q z = {E' r Bg - E' r B g ) /p 0 + (k zr dT e /dr) + {k zz dT e /dz) . 


In the case of an applied magnetic field, the components B r and B z are the sum of the steady-state applied field and an 
induced field that is initially time-varying but eventually converges to a steady-state value. The set of MHD equations 
are closed by the generalized Ohm’s law. 


E' = E + (u x B) = r?j + 


(j x B) - Vp e 

en e 


(7) 


and Ampere’s law. The models for transport coefficients, equation-of-state, ionization and thermal nonequilibrium 
are described in Ref. [26], Implementation of the entire equation set given in Eqns. (6) represents a significant en- 
hancement in simulation capabilities relative to the previous work by Sankaran et al. 24-26 where B r , B z , and u e , and 
non-conserved, secondary (calculated) quantities je and Eg, were all zero. 
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B. Finite Volume Formulation 


The solution of Eqs. (6) in a conservative manner is performed by writing each equation in their integral form. All 
eight conservation laws (mass, three components of momentum, three components of magnetic flux, and total energy) 
can be individually integrated over the cell area. 



( 


dFr dFA 

dr dz ) 


4- Source 


dA. 


A vector F = [F z , — F r ] is defined such that 


VxF|„ = f: 


dr dr dz 


(8) 


(9) 


This definition allows us to write the right hand side of Eq. (8) as the curl of F. Using Stokes’ theorem, the curl over 
the cell area can be expressed as the line integral around the edges of the cell. 



^dA = j F dl + JJ (Source) dA 


( 10 ) 


Since this equation is true irrespective of the shape of the control volume, the fluxes F r and F z can written in the (r, z) 
coordinate system, as shown in Eqs. (6), regardless of the shape of the cell. 


C. Magnetostatic Model 



A separate magnetostatic solver 31 is used to model the 
applied magnetic field topology of the experiment being 
simulated. The magnetic coils in these experiments have 
tens of thousands of loops, while the coils in magneto- 
static simulation had only a handful of loops. Conse- 
quently, the primary goal of our magnetostatic modeling 
effort was to recreate a “reasonable/realistic’’ topology. 

We have modeled the applied magnetic field from the 
magnetic nozzle experiments described in Refs. [5-15, 

22, 32]. In Fig. 1 we show the nozzle topology for the 
experiment described by Kuriki and Okada 6 generated 
using the magnetostatic solver. The computationally- Figurc \ MagnetostaUe simulation or flux lines in the r - , plane of 
, , c , the nozzle in Ref. [6]. 

generated held was compared to published values from 

the experimental setup. Despite our simplifying the 

problem by using only a few magnet coil loops, the magnetostatic solver produced a magnetic field topology pos- 
sessing field values that were within a few percent of those measured. 6 

The magnetostatic solver 31 outputs the relevant variables (B r , B z ) on an orthogonal (r, z) grid. An interpolation 
script was used to transfer the magnetic field data onto the non-orthogonal grid used by the MHD code. 


IV. Code Validation 


Under the assumptions of very low magnetic Reynolds number and isotropic, constant plasma resistivity, Faraday’s 
law reduces to a magnetic diffusion problem governed by the equation 


— = —V 2 
dt 


( 11 ) 


We consider the case of a plasma in the presence of a constant solenoidal field in the r — z plane. At time t = 0 a non- 
uniform current j = j z (r) z is introduced into the plasma. The current resistively diffuses until it is uniform, and the 
induced magnetic field (in the 0 direction) due to the current relaxes to its equilibrium value. Assuming axisymmetry 
and no variation of Bg in the 2 -direction, the resistive diffusion of the azimuthal, induced magnetic field is given by 

dt Ho \ dr 2 r dr r 2 / 


5 of 8 


American Institute of Aeronautics and Astronautics Paper 2007-4386 



We suppose that the initially non-uniform current produces a magnetic field 

= <13 > 

where R is the outer radius of the domain and l is the total current. The boundary conditions on B g at r = 0 and 
r = R are 


B*(0,i) = 
Bg{R,t) = 


0, 

1^0 1 
2ttR' 


The time-dependent solution for the induced magnetic field is 


(14) 
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Pn = {f 0 9{r) - r - Jl & r ) dr ) / (j 0 Jl (t r )' r ‘ Jl (| r ) dr ) ! 


(17) 


and 

g(r) = B e (r,0)-^r. (18) 

Here J\ is the Bessel function of the first kind and z n are its 
zeros. We found that the summation given in Eq. (16) could 
be truncated at n = 5 without a loss in accuracy. 

In our test case, I = 10 4 A, R = 0.05 m and r] = 
10 -4 SI. m. The plasma density was chosen to be high 
(n e = m = 5 x 10 22 m -3 ) so that diffusion processes 
would dominate convection and would thus enable us to 
compare the results of the simulation with diffusion theory. 
We simulated the test case by itself, against a uniform axial 
magnetic field background, and against the solenoidal ap- 
plied magnetic field (with components B r and B z ) given in 
Fig. 1. The solution in Eq. (16) is unaffected by an applied, 
purely axial magnetic field, and the numerical solutions ex- 
hibit good agreement in this case. For the solenoidal field, 
the radial component is exactly zero at the axial plane that 
bisects the ring coil. A comparison between the numerical 
results for Be on the axial plane where B r — 0 and the an- 
alytical solutions at different instances in time are presented 
in Fig. 2. 

From the figure, we observe that the code accurately 
models the time-dependent evolution of the azimuthal mag- 
netic field and the magnetic field components B r and B z 
(from Fig. 1) remain unperturbed. While the temporal evo- 
lution of Be along the axial plane bisecting the ring coil was 
generated against the solenoidal applied field, it is worth 
noting that the numerical and analytical solutions were in 
agreement for the other two applied magnetic field test cases 
(no field and a purely axial applied field). 

Another important result is that the code handled a 
three-dimensional B-field without noticeable deviation in 
the value of V • B from zero. Without any assumptions on 
V • B, Faraday’s law is, 



Radius [cm] 


Figure 2. Comparison between analytical and numerical solu- 
tions for Be in the magnetic diffusion test case. Numerical solu- 
tions were calculated against the magnetic field given in Fig. 1. 
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( 19 ) 


= —V x E - (V ■ B) u. 
o c 

While it is physically true that V ■ B = 0, it is often not true numerically owing to truncation and round-off errors. 
In fact, this has been a source of difficulty in extending MHD codes to three-dimensional vector fields since these 
small deviations in the divergence can lead to numerical instabilities that are unbounded. 33 Consequently, a significant 
amount of work has gone into developing MHD solution algorithms that strictly enforce the condition V • B e 0, In 
our case, because the second term on the right-hand side of Eq, (19) is small, owing to the fact that both V • B and |u| 
are small, truncation errors do not significantly affect the solution to the diffusion problem. 

Many previous MHD codes 34,35 have employed the ideal-MHD assumption because they were developed for the 
purpose of simulating problems with high magnetic Reynolds number - either due to very low resistivity (as is the 
case in fusion studies) or due to the dominance of convection (as is the case in many space and astrophysical plasmas). 
In such cases the absence of any significant dissipation allows for the growth of numerical oscillations 33 due to errors 
in the calculation of V • B. Fortunately, in many plasma thrusters and magnetic nozzles employed for propulsive 
purposes, the magnetic Reynolds number Re m ~ 0(1) and the plasma’s resistivity provides a dissipative mechanism 
within the domain that can damp the numerical instabilities in much the same way that viscosity damps instabilities in 
the solution of the momentum equation. 


V. Conclusions and Future Work 

We have presented the formulation of a simulation tool that is capable of supporting multidimensional numerical 
studies of plasma flows and detachment in magnetic nozzles. The code is based on a conservation formulation of the 
MHD governing equations and contains realistic models of various physical and transport phenomena. The applied 
magnetic field solution is generated using a separate magnetostatic solver. This field is used as an input for the plasma 
flow modeling. There are several proposed mechanisms that could drive the detachment process; the continuum-based 
formulation of the numerical problem limits us to the study of the resistive and kinetic detachment mechanisms, The 
code has demonstrated the ability to model a plasma problem in which all three components of the magnetic field were 
present, and during this simulation deviations from the condition V • B = 0 were small. 

In a future paper, we shall describe simulations of a plasma flow in a magnetic nozzle. For validation purposes, we 
shall base our simulations on experiments described in Refs. [5-15,22,32]. Aside from code validation, the primary 
purposes of these and other future simulations are to: 

• Investigate the feasibility of contouring the magnetic fields such that the plasma accelerates to super-Alfvenic 
speeds and detaches, 

• Vary the plasma properties to investigate whether resistive or kinetic detachment occurs, 

• Generate magnetic field topologies that might more readily lead to plasma detachment, 

• Prescribe guidelines and scaling parameters with relevance to spacecraft propulsion magnetic nozzle design. 
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